%\documentclass[preprint,amsmath,aps,superscriptaddress,showpacs,a4paper]{revtex4}
%SM\documentclass[twocolumn,12pt,cites,graphicx]{article}
\documentclass{article}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{bm}
\usepackage{booktabs}
\usepackage{float}

%\begin{document}

%-----------------------------------------------------------------------

\title{Lattice-Boltzmann simulations of flow in porous media with oilfield applications}

\author{J.\ Yang and E.\ S.\ Boek\\ %\thanks{e-mail:e.boek@imperial.ac.uk} \\
\small Qatar Carbonates and Carbon Storage Research Centre (QCCSRC)\\ 
\small Department of Chemical Engineering\\ 
\small Imperial College London\\ 
\small South Kensington Campus\\ 
\small London SW7 2AZ\\ 
\small United Kingdom\\
\small \texttt{e.boek@imperial.ac.uk}
}

%\author{Dagwood Remifa\thanks{Thanks to the editors of this wonderful journal!}\\
%\small Department of Inconsequential Studies\\[-0.8ex]
%\small Solatido College, North Kentucky, USA\\
%\small \texttt{remifa@dis.solatido.edu}\\
%\and
%Forgotten Second Author\\
%\small School of Hard Knocks\\[-0.8ex]
%\small University of Western Nowhere\\[-0.8ex]
%\small Nowhere Uvherdov\\
%\small \texttt{no1remembers@me.woe.edu}
%}

%\date{\dateline{Jan 1, 2009}{Jan 2, 2009}{Jan 3, 2009}\\


\begin{document}
\maketitle

\begin{abstract}
In this paper, we present a new computational fluid dynamic (CFD) method: lattice Boltzmann method (LBM) to predict the transport properites of porous media. We describes the details on the LBM method, including two sinlge phase LB model using Single-Relaxation-Time (SRT) scheme and Multi-Relaxation Time (MRT) sheme, three multi-component LB models for binary flow simulations. Several numerical experiments are carried out to veryfy the implementation. A permeability calculation for porous medium with cylinders with elliptical cross section are carried out, simulated premeability agrees well with analytical solutions. The multi-component LB models are used to study several fundermental problems for complex flow in porous media including phase separation, laplace law, effect of wetting surface and poisseuille flow for a binary fluid with viscosity contrast, we demonstrate that the model can accurately reproduce these phenomenas. 

    The lattice Boltzmann method is applied to study two improtant mutli-phase phenomenas: capillary fingering and snap-off. As enhanced oil recovery (EOR) normally invlves displacement of one phase by injection of another fluids phase, stability of this fluid-fluid displacement affects the degree of EOR (capillary fingering), at the same time, a big quantity of oil remains in an oil reservoir after waterflooding due to the competation of capillary between phases (snap off). The capillary fingering and snap-off might lead to inefficient recoveries, therefore it is importnat to understand the mechanisms.The lattice Boltzmann method is used to study these phenomenas, and successfully replicate capillary fingering snap-off phenomenas, good agreements with literature data are obtained. And we correctly observe that the arise of snap-off phenomena depend on the strength of surface tension. It is concluded that the lattie Boltzmann method is a promising tool for predicting the transport properties of porous medium.    
\end{abstract}

%SM:
%\author{J.\ Yang and E.\ S.\ Boek
%%\thanks{Present address: Insert New Address Here.}
%\\
%Qatar Carbonates and Carbon Storage Research Centre (QCCSRC),
%Department of Chemical Engineering,
%Imperial College London,
%South Kensington Campus,
%London SW7 2AZ
%United Kingdom
%Address, Town, Country. 
%E-mail: e.boek@imperial.ac.uk 
%}

%\date{Received XXXXth Month, 200X\\Accepted XXXXth Month,
%200X\\DOI: 10.1039/}

%\begin{document}

%\author{J.\ Yang}
%\affiliation{
%Qatar Carbonates and Carbon Storage Research Centre (QCCSRC), 
%Department of Chemical Engineering, 
%Imperial College London, 
%South Kensington Campus, 
%London SW7 2AZ
%United Kingdom
%}
%
%\author{E.\ S.\ Boek}
%\affiliation{
%atar Carbonates and Carbon Storage Research Centre (QCCSRC),
%Department of Chemical Engineering,
%Imperial College London,
%South Kensington Campus,
%London SW7 2AZ
%United Kingdom
%}
%\email{e.boek@imperial.ac.uk}

\date{\today}
%-----------------------------------------------------------------------

%\begin{abstract}
%Abstract
%\end{abstract}

%\pacs{}

%\maketitle

%\begin{document}

%\section{Abstract}

\section{Introduction}
Much research in recent years has focused on the prediction of transport properties of fluids in reservoir rocks including absolute permeability, relative permeability and capillary pressure using numerical methods. This is because numerical methods can predict the transport properties of reservoir rocks with higher efficiency (in several days) and lower cost%(less than \$10000)
compared with experiments which might take several months and therefore would cost significantly more.%and up to \$100000.
The lattice Boltzmann method (LBM)\cite{doolen_1990}\cite{chen_1992} is a novel Computational Fluid Dynamic (CFD) technique which is able to handle the flow through extremely complex geometry without simplification \cite{cancelliere_1990}. Whereas traditional CFD techniques solve the Navier-Stokes equations directly, based upon the conservation of the macroscopic properties in a fluid, the LBM models fluid flow based on the evolution of microscopic Particle Distribution Functions (PDFs). %the distribution function of discrete pseudo-particles.
It has been shown that the Navier-Stokes equations emerge and are recovered from the LBM equations at sufficiently large length scales.~\cite{succi_book}
The motion of these PDFs is confined to the nodes of a specified grid, making it possible to simulate the behaviour of single and multi-phase fluids in and around complex geometries\cite{boek_2010}. This method is easier to implement than other conventional CFD technique and is ideal for parallel computing due to the fact that most of the operations in LBM are local. The LBM has many possible uses in the petroleum physics field and has been used to study the transport properties of porous media \cite{jin_2004}\cite{Pan_2006}\cite{pan_2004}\cite{ramstad_2009}\cite{boek_2010}.

\section{Methodology}
\subsection{Single Relaxation Time BGK(Bhatnagar-Gross-Krook) lattice Boltzmann Method}

%\subsection{Single-phase lattice Boltzmann method}
The lattice Boltzmann Method (LBM) is a special discretization of the Boltzmann-BGK equation. Discretization of space, velocity and time are carried out in LBM. This procedure greatly simplifies the original Boltzmann equation. The location of particle distribution functions (PDFs) in space is restrained to the nodes of a lattice grid, and the pdf velocities are simplified to a very limited number of lattice velocities. We take a 2-D model as an example. This model is well known and widely used in the application. It is two dimensional and contains 9 velocities with the name D2Q9. This model has been proposed by Qian et al\cite{qian1992}. In LBM, we assume that all the PDFs have the same uniform mass (normally 1 mass unit is taken for simplicity). The lattice unit($lu$) and time steps ($ts$) are important length and time units in LBM. We only discuss uniform mesh in this chapter ($\Delta x= \Delta y$).


\begin{figure}[H]
\begin{center}
\includegraphics[width=2in]{d2q9.eps}
\end{center}
\caption{D2Q9 lattice and velocities}\label{d2q9}
\end{figure}

Figure \ref{d2q9} shows the discretized velocity space $\{ \mathbf{e}_i \} ,\quad (i=0..8)$. The lattice velocity can be written as:

\begin{equation}
\mathbf{e}=e\begin{bmatrix}
0 & 1 & 0 &-1 &0 &1& -1& -1& 1 \\
0 &0 &1 &0 &-1& 1& 1& -1& -1
\end{bmatrix}
\end{equation}

where $e=\Delta x / \Delta y$ is the local lattice speed, and has is related to sound speed as $c_s=\frac{e}{\sqrt{3}}$

Correspondingly, the continuous distribution functions associated with velocity are written as $f_i(x,t),\quad (i=0..8)$. We can obtain the Lattice Boltzmann equations for the D2Q9 model (Single Relaxation Time BGK) as:

\begin{equation}
f_i(\mathbf{x}+\mathbf{e}_i \Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{f_i(\mathbf{x},t)-f^{eq}_i(\mathbf{x},t)}{\tau}
\end{equation}

Collision of the particles can be considered as a relaxation process towards equilibrium. The equilibrium distribution function is defined as:

\begin{equation}
f^{eq}_i(x)=w_i\rho (x)[1+3\frac{\mathbf{e}_i \cdot \mathbf{u}}{e^2}+\frac{9(\mathbf{e}_i \cdot \mathbf{u})^2}{2e^4}-\frac{3\mathbf{u}^2}{c^2}]
\end{equation}

where the weight coefficients for the D2Q9 model are:
\begin{equation}
w_i=\begin{cases}
4/9 & i=0 \\
1/9 & i=1..4 \\
1/36 & i=5..8
\end{cases}
\end{equation}

The macroscopic transport equations for mass, momentum and energy can be derived from the Boltzmann equation using a Chapman-Enskog expansion.\cite{ekexpansion} The kinematic viscosity $\nu$ in the D2Q9 model is obtained as:
\begin{equation}
\nu=c^2_s(\tau-\frac{1}{2})\Delta t
\end{equation}

%Its unit are $lu^2ts^{-1}$.
Note that $\tau>1/2$  for positive viscosity. Numerical difficulties can arise as $\tau$ approaches $1/2$. A value of $\tau=1$ is safest and leads to $\nu=1/6 lu^2ts^{-1}$.\cite{LBMODELING} The pressure is given by the equation of state for an ideal gas:
\begin{equation}
p=\rho c^2_s
\end{equation}


To implement a Lattice Boltzmann simulation code, four major steps should be included:
\begin{itemize}
\item Initialisation of distribution function $f_i(\mathbf{x},0)$
\item Collision step
	\begin{equation}
	f'_i(\mathbf{x},t)=f_i(\mathbf{x},t)+\Omega_i(\mathbf{x},t)
	\end{equation}
\item Streaming step
	\begin{equation}
	f_i(\mathbf{x}+\mathbf{c}_i\Delta t,t+\Delta t)=f'_i(\mathbf{x},t)
	\end{equation}
\item Computation of Macroscopic hydrodynamic quantities
	\begin{eqnarray}
	&\rho(\mathbf{x},t)=\sum_{i} f_i(\mathbf{x},t)& \\
	&\rho \mathbf{u}(\mathbf{x},t)=\sum_{i} \mathbf{c}_i f_i(\mathbf{x},t)&
	\end{eqnarray}
\end{itemize}

\section{Bounceback Boundary Conditions}

The bounceback condition plays a major role in the LBM simulation due to its simplicity, versatility and powerful capability of dealing with extremely complex boundaries. This boundary condition is usually used at fluid-solid interfaces due to its correspondence to the no-slip condition. This boundary condition is illustrated in Figure \ref{bounceback}. The densities moving toward the solid are bounced back into the fluid domain with the incoming directions. In the D2Q9 model, the bounceback condition can be described in terms of equations as:

\begin{eqnarray}
f_2=f_4,&f_5=f_7,&f_6=f_8
\end{eqnarray}

The standard bounceback condition places the boundary on the lattice nodes. The bounceback scheme conserves mass and momentum, but the accuracy is first order, while LBM is of second order. Inamuro \cite{inamuro1995} found that the error produced by the bounceback condition is sufficiently small if the relaxation parameter $\tau$ is close enough to 2. The bounce back conditions can be used without any influence on the order of the LBM, if $\tau$ is chosen within a suitable range. Furthermore, the bounce back condition is the most efficient one for arbitrarily complex geometries (Breuer et al. 2000)\cite{breuer2000}. Many researchers have contributed to this ongoing discussion \cite{noble_1995}\cite{ziegler1993}. A second order scheme with the name Half-Way bounce back condition was proposed by Ziegler\cite{ziegler1993}. In this boundary condition, the surface is a solid boundary placed between two neighbouring lattice sites with the same distance $\Delta x /2$. It is illustrated in Figure(\ref{bounceback}).


\begin{figure}[H]
\begin{center}
\input{bounceback.pstex_t}
\end{center}
\caption{Bounceback condition}\label{bounceback}
\end{figure}

\section{Multi-Relaxation-Time (MRT) scheme for the lattice Boltzmann method}
To overcome the limitations of the single relaxation time LBGK model, the Multi-Relaxation-Time (MRT) scheme was introduced. This scheme allows independent adjustment of bulk and shear viscosities which significantly improves the numerical stability for a low viscosity fluid. In the LBGK model, the collision term is relaxed by single parameter $\tau$, while in MRT it is relaxed by a matrix $\Lambda$  \cite{dhumieres1992}\cite{dhumieres2001}:

\begin{equation}
f_i(x+c_i \Delta t)-f_i(x,t)=-\Lambda_{ij}[f_j-f_j^{eq}], \quad i=1,2,\ldots,b
\end{equation}

The matrix $\Lambda$ is a full matrix with constant ?. The LBGK model can be obtained by specifying $\Lambda$ as a diagonal matrix with identical values:

\begin{equation}
\Lambda_{ij}=\frac{1}{\tau}\delta_{ij}
\end{equation}

The macroscopic quantities are calculated in the same way as in the LBGK model. Instead of considering distribution functions, MRT employs several moments corresponding to macroscopic quantities and their flux. These quantities can be relaxed with different time scales. A matrix $M$ transforms the distribution functions $f_i$ from the distribution space to the moment space\cite{dhumieres1992}\cite{dhumieres2001}\cite{lallemand2000}:

\begin{equation}
m=M\cdot f, \quad f=M^{-1}\cdot m
\end{equation}

The collision is carried out in the moment space by multiplying the transformation matrix $M$; the left and right hand side of equation \ref{collisionmrt} can be transformed into the moment space as:

\begin{equation}
f_i'(x,t)=f_i(x,t)-\Lambda_{ij}[f_j-f_j^{eq}]
\label{collisionmrt}
\end{equation}
\begin{equation}
\mathbf{m}'=\mathbf{m}-\mathbf{S}[\mathbf{m}-\mathbf{m}^{eq}]
\end{equation}

where $\mathbf{m}^{eq}=\mathbf{M}\mathbf{f}^{eq}$ is the equilibrium equation in moment space. $S=M \Lambda M^{-1}=diag(s_1,s_2,\ldots,s_b)$. The corresponding relaxation time for moment $m_i$ is $s^{-1}_i$. After the collision step, the moment $m'$ is transformed back into distribution function space by multiplying by $M^{-1}$ for the streaming step which will be carried out in the same way as in the  LBGK model. The moment space for the D2Q9 model is:

\begin{equation}
m=(\rho,e,e^2,j_x,q_x,j_y,q_y,p_{xx},p_{xy})^{T}
\end{equation}

where $e$ is the energy, $j_x,j_y$ are the components of the momentum in $x$ and $y$ directions, $q_x,q_y$ are energy fluxes, $p_{xx},p_{xy}$ is the stress tensor. The relaxation parameters and equilibrium functions of the moments are:



\begin{equation}
S=(0,s_e,s_{e^2},0,s_q,0,s_q,s_{\nu},s_{\nu})
\end{equation}
\begin{equation}
m^{eq}=\rho(1,-2+3u^2,\alpha+\beta u^2,u_x,-u_x,u_y,-u_y,u^2_x-u^2_y,u_xu_y)^{T}
\end{equation}

where $\alpha$ and $\beta$ are parameters, and kinematic viscosity and volume viscosity are:
\begin{eqnarray}
&\nu=c^2_s(\frac{1}{s_{\nu}}-\frac{1}{2})\Delta_t& \\
&\zeta=c^2_s(\frac{1}{s_{e}}-\frac{1}{2})\Delta_t&
\end{eqnarray}

\section{Multi-component lattice Boltzmann method}
A multi-component system consists of separate
chemical components such as oil and water. Due to their economic
importance, such systems have been studied extensively. Three different multi-component lattice Boltzmann models are implemented, and their performance is evaluated with some numerical experiments.

\subsection{Shan-Chen pseudo potential model}
A pseudo potential model for multi-phase and multi-component lattice Boltzmann simulation was proposed by Shan and Chen\cite{shanchen}. The principal characteristic of this model is an interaction force between particles which is introduced to have a consistent treatment of the equation of state for a non-ideal gas. For the D2Q9 model,the attractive force between the PDFs is described as:

\begin{equation}
F(\mathbf{x},t)=-G\Psi (\mathbf{x},t)\sum_{i=1}^{8}w_i \Psi (\mathbf{x}+\mathbf{e}_i\Delta t,t)\mathbf{e}_i
\label{Shanchenforce}
\end{equation}

where $G$ is a parameter determining the interaction strength between neighbouring particles. It also determines whether the interaction is attractive or repulsive. In multi-component flow, same component particles will attract each other while particles from different components will repel each other. From this, phase separation can be obtained. And $\Psi$ is the interaction potential with the form:

\begin{equation}
\Psi (\rho)=\Psi_0exp(\frac{-\rho_0}{\rho})
\label{potentialfunction}
\end{equation}

where $\Psi_0$ and $\rho_0$ are constant parameters.

\subsection{The Free Energy lattice Boltzmann model}
Swift et al \cite{freeenergy} developed a new thermodynamically consistent
binary fluid LB model by introducing an equilibrium state associated
with a free energy functional, corresponding pressure tensors and
chemical potentials. A correct choice of collision rules ensures that
the system evolves toward minimization of the free energy
functional.

The thermodynamic properties of a binary fluid system can be
described by a Laudau free energy functional:

\begin{equation}
\Psi=\int_V (\psi_b+\frac{\kappa}{2}(\partial_{\alpha} \phi)^2)dv
+\int_S \psi_s ds
\end{equation}

where $\psi_b$ is the bulk free energy density and has a form:
\begin{equation}
\psi_b=\frac{c^2}{3} \rho ln\rho +A (-\frac{1}{2}\phi^2+
\frac{1}{4}\phi^4)
\end{equation}

$\phi$ is the order parameter representing the concentration of
components defined as:

\begin{equation}
\phi=\frac{\rho_a-\rho_b}{\rho_a+\rho_b}
\end{equation}

The hydrodynamics and thermodynamics of the binary fluids are
described by the Navier-Stokes equations and a convection-diffusion
equation:

\begin{equation}
\frac{\partial \phi}{\partial t}+\nabla (\phi u)= \mu \nabla ^2 \mu
\end{equation}

A new distribution function $g_i(x,t)$ is introduced to describe the
concentration $\phi=\sum\limits_i g_i$ and is related to the
convection and diffusion. The distribution function $f_i(x,t)$ is
related to the fluid density and momentum as usual.

The time evolution equations using the MRT scheme can be described as:
\begin{eqnarray}
\text{Collision Step:} & f'_i(x,t)=f_i(x,t)-M^{-1}SM(f_i-f^{eq}_i) & \\
& g'_i(x,t) =g _i^{eq}(x,t)& \\
\text{Streaming Step:} & f_i (x+e_i \Delta t,t+\Delta t) =f'_i(x,t)&
\\
&g_i (x+e_i \Delta t,t+\Delta t) =g'_i(x,t)&
\end{eqnarray}

An appropriate choice of the equilibrium distribution functions can
reproduce the macroscopic equations in the continuum limit \cite{kusumaatmaja_phd}.

\subsection{The Colour Gradient lattice Boltzmann Model}
An immiscible fluid model developed from Lattice Gas Cellular
Automata was introduced by Gunstensen and Rothman \cite{color1991}. The
particles in this model are coloured either red or blue and therefore it
is normally called the Colour Gradient Method. The surface tension is
introduced by adding a perturbation to the collision operator \cite{color1991} while
adhering to the Navier-Stokes equations in homogeneous
regions. A recolouring step is involved after the surface tension
perturbation calculation in order to achieve zero diffusivity of one
colour into the other.\\

We use $f_i^r$, $f_i^b$ and $f_i$ to denote the distribution functions
of a red fluid, a blue fluid and their combination respectively. A perturbation is computed to generate the surface tension. The
surface tension can be expressed as a local anisotropy in the
pressure: the pressure measured normal to the surface is larger than
that tangential to the surface. The pressure in a LBM simulation is proportional to the density ($p=c_s^2\rho$), so surface
tension can be generated by preferentially placing particles in
directions normal to the interface rather than tangential. Mass and
momentum should be conserved. The color gradient $G$ is defined as:

\begin{equation}
G(x,t)=\sum\limits_i e_i \sum\limits_j(f_j^r(x+e_i\Delta
t,t)-f_j^b(x+e_i\Delta t,t))
\end{equation}

The perturbation of the populations is:
\begin{equation}
f''_i(x,t) = f'_i(x,t) +A|G(x,t)|cos(2\theta)
\end{equation}


A redistribution of colour forces the particles to move towards the
regions occupied by particles with same colour. This recolouring step
enables us to achieve the separation of the fluids. It is carried out
by the following maximization problem \cite{color1991}:

\begin{equation}
W(f_i^{r''},f_i^{b''})=\max\limits_{f_i^{r''},f_i^{b''}}[\sum\limits_i
(f''_{r_i}-f''_{b_i})e_i]
\end{equation},

subject to the following constraints:

\begin{eqnarray}
&\rho''_r=\sum\limits_i f_i^{r''}=\rho_r&\\
&f_i^{r''}+f_i^{b''}=f_i&
\end{eqnarray}

A general colour gradient lattice Boltzmann model can be summarised
as:

\begin{itemize}
\item Single phase collision.
    %\begin{equation}
    %f_i=f_i^r+f_i^b
    %\end{equation}
\item Add a surface tension perturbation to $f'_i$ obtaining $f''_i$
\item Recolouring
\item Streaming
\end{itemize}

\section{Verifications}
\subsection{2D Flow Past a Cylindrical Body}

\subsection{Simulation of flow in fibrous porous media}

A porous medium with rectangular periodic arrays of cylinders with elliptical cross section \cite{yang_2000} is investigated to evaluate the performance of our code. A square packing configuration for elliptical rows of fibers is shown in Figure\ref{fibrous_geo}. Periodic boundary conditions are applied in the $x$, $y$ and $z$ directions. An analytical formulation of the permeability for this geometry is derived by Phelan et al.(1996)\cite{phelan_1996}, which is used to evaluate the simulated permeability. The analytical expression is as follows: \\

\begin{equation}
\end{equation}

\begin{figure}[H]
\begin{center}
\includegraphics[width=1.0in]{fibrous_porous_geo.eps}
\end{center}
\caption{Geometry of fibrous porous media}\label{fibrous_geo}
\end{figure}

The permeability is calculated using a 150x150x10 mesh in the simulation. The results are compared with both LB results from a previous study \cite{Keehm_PHD}, using a Single Relaxation Time (SRT) alorithm, and the analytical solution. It can be seen from Figure (\ref{elliptical_per}) that good agreement with the analytical solution is obtained. In fact, our solution seems to approximate the analytical solution more closely than the previous LB results. This is probably due to the the fact that we have used a Multiple Relaxation Time (MRT) scheme, which gives more accurate results than the SRT scheme~\cite{Keehm_PHD}. In low porosity porous media simulations, the prediction of the permeability from LBM showed a bigger error \ref{eclliptical_per_error}. This is due to the lack of fine mesh between adjacent cylinders.

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{cylinder_permeability.eps}
\end{center}
\caption{Simulated dimensionless permeability of arrays of cylinders with elliptical cross sections. The red solid line shows analyticao solution, blud crossing points are simulated permeability from our MRT-LBM, the black trangle points are literature data from Keehm's paper\cite{Keehm_PHD}}.  \label{elliptical_per}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{cylinder_per_error.eps}
\end{center}
\caption{Relative error of simulated permeability.The blue points line shows error of our MRT-LBM, the black points line shows the error of Keehm's results. Our results show a decreasing relative error with the increase of porosity, while Keehm's results yield big errors on both low and high porosity porous medium}. \label{eclliptical_per_error}
\end{figure}
\subsection{Phase separation}
 Using the Shan-Chen pseudo potential model, we can simulate liquid-vapour phase separation and its dynamics. A $200 \times 200$ domain with average density of $200 mu lu^{-1}$ is computed. The initial density is specified with random variations to prevent a metastable situation. We used the parameters $\Phi _0=4$, $\rho_0=200$ and $G=-120$ for the calculation of interaction force term. The surface area is minimized as a consequence of free energy minimization. In the liquid-vapor system, condensation and evaporation phenomena are also discovered during the simulation. This may be due to the relatively high vapor density.

\begin{figure}[H]
\begin{center}
\includegraphics[width=5in]{phaseseperation.eps}
\end{center}
\caption{Liquid-vapour phase separation dynamics at times 0, 100, 200,400,2000,6000 $ts$}
\end{figure}

\subsection{Maxwell Construction}
Flat interfaces are exceptionally important because the vapour pressure above the interface at equilibrium is the saturation vapour pressure. The pressure across the flat interface should be zero according to the Laplace law:\begin{equation}\Delta p=-\gamma(\frac{1}{R_1}+\frac{1}{R_2})\label{lap_law}\end{equation} Where $\gamma$ is the surface tension, $R_1$ and $R_2$ are principle radii curvatures. The density of the liquid and the vapour from the simulation can be compared with the Maxwell Construction result \cite{LBMODELING}. A 50*200 domain was set to compute the yield density from Shan-Chen model. A potential function in equation (\ref{potentialfunction}) with $\Psi_0=4, \quad \rho_0=200$ and $G=-120$ is computed. The Shan-Chen multiphase model is used to separate the phases. We use two different force terms: the first is the shifting velocity scheme proposed by Shan and Chen \cite{shanchen1995}, which has been widely used in LBM multiphase simulations; another force term is presented by Guo et al.\cite{guo2002} with consideration of discrete lattice effects. The simulation results are compared with the numerical solution from the Maxwell Construction and given in the Table below.

\begin{table}[!htbp]
\begin{center}
\begin{tabular}{lccc}
\toprule

& $\rho_{liquid}$ & $\rho_{vapor}$ & $|u_{max}|$  \\
\midrule
Shan-Chen & 524.4 & 85.75 & 1.6291E-4 \\
Guo force term & 514.7 &79.69 & 1.1204E-4 \\
Analytical solution\cite{LBMODELING} &514.64 & 79.705 & \\
\bottomrule

\end{tabular}
\end{center}
\caption{Liquid and vapor density yield from Shan-Chen model with two different force terms and the numerical solution from the Maxwell Construction. All values are in lattice units}
\end{table}

The results show that liquid and vapour density yield from the Shan-Chen model with the original force term has a significent difference from the analytical solution, while the Guo force term leads to very accurate result. The way of incorporating the force term also has big effects on the spurious velocity $u_max$. The Guo force term also gives a smaller spurious velocity across the interface. Guo et al. \cite{guo2002} showed that the Shan-Chen original force term is a first order scheme, which can only recover the mass equation, and has an error term $\nabla \cdot[\rho \nu^2(\nabla a+\nabla a^T)]$ in the momentum equation. Guo's force term, on the other hand, can recover both the mass and momentum equations. This may be the reason leading to a more accurate liquid and vapour density compared with the Shan-Chen force term.


\subsection{Surface tension}
By simulating a series of drops of different size and measuring the inside and outside pressure, the surface tension can be estimated with the Laplace law, see equation (\ref{lap_law}). The slope of 1/radius vs. pressure difference will be the surface tension. The following numerical examples use an improved Shan-Chen model (\ref{expansionforce1a}) which enables the surface tension to be changed independently of the EOS.

\begin{equation}
F=-c^2_sG\Psi \nabla \Psi -\frac{1}{2}\Delta x^2c^2_sG\Psi \nabla \nabla^2 \Psi +o(\Delta x^3)
\label{expansionforce1a}
\end{equation}

A non-ideal EOS(\ref{sukopeos2}) with the potential function defined in Equation(\ref{potentialfunction}) is used.
\begin{equation}
p(\rho)=\rho c^2_s +\frac{1}{2}Gc^2_s \Psi ^2(\rho)
\label{sukopeos2}
\end{equation}

We now present the pressure fields computed for different sizes of drops and different values of $\kappa$.

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{surface_tension.eps}
\end{center}
\caption{Surface tension estimation with different $\kappa$}
\label{sur1}
\end{figure}


\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{surface_tension2.eps}
\end{center}
\caption{Surface tension $\sigma$ change along with parameter $\kappa$}
\label{sur2}
\end{figure}

Fig (\ref{sur1}) shows that the Laplace law is satisfied with $\kappa=12,20,30$. Different values for $\kappa$ can result in different values of the surface tension $\sigma$, calculated using the Laplace law. In general, we observe that results for small surface tension fit the Laplace law well, while a larger surface tension produces a bigger error. Fig (\ref{sur2}) shows the linear relation between the parameter $\kappa$ and the surface tension $\sigma$. This numerical example demonstrates the validation of the new Shan-Chen model with improvements which enable us to control the surface tension seperately from the equation of state.

\subsection{Simulation of wetting surface}To verify that the lattice Boltzmann method is able to model different wetting surfaces, simulations for droplets on a wetting surfaces with different contact angles were performed. The red fluid was initialized as a cubic configuration, then the code was performed until the system reach equilibrium. Droplets on a wetting surface with contact angles $20^\circ,90^\circ,150^\circ$ are simulated with the Free Energy model and are shown in Figure(\ref{wettingb}).

 \begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{wetting.eps}
\end{center}
\caption{Droplets simulation on wetting surface with contact angle $20^\circ$(left),$90^\circ$(middle) ,$150^\circ$ (right)}\label{wettingb}
\end{figure}The conatct angle determins the wetting perferences of the solid, which is critical to EOR processes. Flooding technology is used widely in EOR processes, the degree of wettability of the solid affect the injection of water into a porous matrix and so the amount of oil in place that will be displaced by waterflooding. The lattice Boltzmann method provide a reliable and efficient way to model wetting surfaces.


\subsection{Poiseuille flow simulation for a immiscible binary fluid system with viscosity contrast}
Simulations of two fluids with a viscosity ratio of 100 in a channel
have been carried out and are compared with theoretical predictions Figure(\ref{CGFE_v}). The densities of the two components are kept the same. The system is initialized with
substance 0 in the middle and substance 1 on both sides near the
walls. An initial density value of 1 is applied to substance 0 and 1. \\

\begin{figure}[H]
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{CG_poisseuille.eps}
%\caption{Simulated velocity of Color Gradient Model \label{CG_v}}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{FE_poisseuille.eps}
%\caption{Simulated velocity of Free Energy Model \label{FE_v}}
\end{minipage}
\caption{Simulated velocity of the Color Gradient Model (left) and the Free Energy Model (right) with viscosity ratio 100} \label{CGFE_v}
\end{figure}

The simulation results of both the Free Energy Model and the Colour Gradient Model
give excellent agreement with the analytical solution.
We should note that the interface thickness of the Free Energy Model is around $6$ lattice
units and the interface location is found at $x=38$ instead of the
initial position of $x=35$. An interface movement of 3 lattice units is discovered from the results. The interface thickness of the Color Gradient Model is around 4, which is smaller than that of the Free Energy Model.
An interface shifting is observed in the colour gradient model of
around 1.5 lattice units which is significantly smaller than the
free energy model. The recolouring algorithm separates the two immiscible fluids very well
and gives nearly 0 diffusivity as expected.


\subsection{Capillary fingering simulation}
Capillary fingering is a kind of hydrodynamic instability which exists in various displacements during oil/gas production. One fluid is pushed by another one, with different viscosity, along a
channel with non-slip walls. This displacement processing can be characterized by a dimensionless number: Capillary number $Ca$. It represent the relative effect of viscous forces versus surface tension acting across an interface between two immiscible liquid. A growing finger of the driving fluid will be produced if the capillary
number $Ca$ is low enough.


\begin{equation}
Ca=\frac{u_t \rho \nu_2}{\sigma}
\end{equation}

where $u_t$ is the velocity of the tip of the finger, $\nu_2$ is the
viscosity of driving fluid and $\sigma$ is the surface tension.


Capillary fingering is investigated by two lattice Boltzmann binary models: the Free Energy Model and the Color Gradient Model. Figure (\ref{fingering2}) shows the evolution of fingers simulated by the Free Energy Model. From top to bottom, they are finger evolution for surface tensions of 0.0678,0.0389,0.01985,0.00992 respectively. No finger will be produced if the surface tension is high. When the surface tension decreases, fingers are observed. However, the smaller the surface tension, the less stable were the interfaces produced. Halpern and Gaver \cite{halpern_1994} studied the fingering phenomenon in a channel.  The width of the fingers produced by different capillary numbers $Ca$ is measured. Their results are plotted in Figure (\ref{fingering1}) with a solid line. A relative width value with the width of the channel (width of the fingering/witdth of the )was used (explain exactly what you mean here). Our results from the Free Energy Model and the Color Gradient Model are shown in the same figure with triangle points and star points respectively. \\

\begin{figure}[H]
\begin{center}
\includegraphics[width=1.5in]{Capillary_Fingering1.eps}
\end{center}
\caption{Fingering evolution with different surface tension at a time interval of 1000 time steps}\label{fingering2}
\end{figure}



\begin{figure}[H]
\begin{center}
\includegraphics[width=2.5in]{fingering.eps}
\end{center}
\caption{Fingering with as a function of Capillary number, the results from the Free Energy model are represented with triangle points, the results from the Color Gradient Model are represented with star points}\label{fingering1}
\end{figure}

Good agreement is achieved although some small discrepancies are found. These discrepancies might come from the boundary conditions. These numerical examples show that both the Free Energy Model and the Color Gradient Model are capable of simulating the capillary fingering phenomenon. It is worth mentioning that, the maximum viscosity ratio for dynamic interfaces simulations is around 20. Knowing this value can help to identify whether the Free Energy Model or the Color Gradient Model is suitable for the specific case. 



\subsection{Simulation of Snap-off Phenomena}

Snap-off is a type of mechanism for the fluid displacement in porous media. \cite{zhao_2010} In a water wetting porous medium which is comprised of pores and throats, the water will accumulate in some regions near the solid boundaries until the oil does not have contact with the solid. Water then fills the throats and separates the oil into droplets. It takes place in the oil migration, and is of great industrial value. LBM has a strong capability of dealing with multi-component flows and therefore we studied this important mechanism in pore scale fluid flow.\\


\begin{figure}[H]
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{snap-off.eps}
\caption{A pore/Throat geometry used in the simulation}\label{snapoffgeo}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{snap-off_2d.eps}
\caption{A pore/Throat geometry in 2D}\label{snapoffgeo_2d}
\end{minipage}
\end{figure}


A geometry shown in Figure(\ref{snapoffgeo}) is established to simulate the snap-off phenomena at the pore scale. It is a tube comprising of equilateral triangles with different inner radii. The solid is set to water wetting  with a contact angle of $60^\circ$. A two dimensional view of the geometry is given in Figure (\ref{snapoffgeo_2d}). Due to the wetting properties, the corner will be occupied by water. This is the region where water accumulates. A binary fluids system which comprises of 1) water of density $1.0$ and 2) oil of density $1.0$ and viscosity of $0.1$ and $0.01$ respectively is simulated. We use the Free Energy Model with a grid 240x40x40, surface tension parameter $\kappa=0.04$, contact angle $60^\circ$ and a body force in $x$ direction of $10^-5$. Initially, the water and oil occupy half of the tube respectively. A body force is applied to simulate the forced imbibition process in the porous medium. The oil will be displaced by the water. Figure(\ref{snapoffs_2},\ref{snapoffs}) shows this displacement process. Water is observed to accumulate in the corners, and invade into the pores. \\

As the oil was displaced by the water, the radii of curvature decreased. The surface tension, determined by the radii of curvature, decreased. As the radii of curvature decreases, the capillary pressure increases. As we can observe in Figure(\ref{snapoffs_2}), the radii of curvature kept decreasing while water was invading the throats.   The contact between volumes of oil in different pores is broken if the surface tension is high enough ($\kappa>0.3$). Some oil will stay in the pore, and will not be displaced by the water.  Figure (\ref{snapoffs_2}) shows clearly this snap off phenomenon. Water is represented by the red fluid, and oil is represented by the blue fluid. Distribution of oil and water is plotted in an interval of 1000 computational time steps. A part of the boundary is removed to better illustrate this processing. \\

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{Snapoff_2.eps}
\end{center}
\caption{Snap off simulation with high surface tension}\label{snapoffs_2}
\end{figure}


Surface tension is the critical condition for the emergence of snap-off. It should reach a critical value, then the oil will be cut into two parts. From the Laplace Law, we can find that the surface tension is controlled by radii of curvature. It is determined by the contact angle and geometry of the pore structure. A simulation with surface tension $\kappa=0.02$ was carried out to investigate the effect of surface tension on snap-off. The simulation results are shown in Figure (\ref{snapoffs}). Snap off phenomena are not found in this simulation, as the capillary pressure determined by the surface tension and radii of curvature of the interface are not strong enough to break the interface in the throat.

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{Snapoff_s.eps}
\end{center}
\caption{Snap off simulation with low surface tension}\label{snapoffs}
\end{figure}



\section{Conclusions}
The lattice Boltzmann method for single phase flow and multi-phase flow simulation offers us an alternative way of predicting the transport properties of flow within porous media. In this study, basic algorithms for single phase three multi-phase LBM models including the Shan-Chen pseudo potential model, the Oxford Free Energy model and the Colour Gradient model were reviewed; to check the validity of the method, we investigated a number of applications, including 1) the permeability of porous media comprised of cylindrical bodies, 2) phase separation, 3) droplets on solid surfaces with different wettabilities and 4) surface tension of a single bubble. The numerical results were compared with analytical solutions and good agreement was obtained. We studied two important phenomena in the displacement of fluids in a reservoir, including capillary fingering and snap-off. For this purpose, we used the Oxford Free Energy LBM. Capillary fingering is an important way of displacing oil in the reservoir, and the role of snap-off is of great importance with regard to the efficiency of oil recovery. Capillary fingering can be observed when a fluid is pushed by a less viscous fluid; good agreement was obtained regarding the finger widths with different capillary numbers with literature data. The snap-off was also successfully captured using the Free Energy LBM. The effect of surface tension on snap-off was also studied; the results indicated that snap-off happens when the surface tension is high enough. A quantitative analysis with theoretical predictions will be future work.\\

These preliminary results provide evidence of suitability of the LBM for flow in porous media and phase behaviour simulation. In contrast to other CFD techniques including the Finite Elements Method and the Finite Volume Methods, this novel CFD technique is of higher efficiency. The inherent parallelism makes the computation on high performance computers easy.\\

We would like to acknowledge that this study is financed by Qatar Carbonates and Carbon Storage Research Centre (QCCSRC). The QCCSRC is funded jointly by Qatar Petroleum, Shell and the Qatar Science \& Technology Park. We thank the High Performance Computing Service of Imperial College London for providing the computing time and technical support.
% for this study.

\bibliographystyle{ieeetr}
\bibliography{ref}

\end{document}
